self-similar propagation of a finite central crack in a finite panel [analogous to the well-known Broberg's problem] (b) stress-wave loading of a stationary central crack in a finite panel

-In this paper, several formulations of moving-singularity finite element procedures for fast fracture analysis are evaluated as to their accuracy and efficiency. INTRODUCTION IN REFS. [1-3] the authors have presented a "moving singular-element" procedure for the dynamic analysis of fast crack-propagation in finite bodies. In this procedure, a singularelement, within which a large number of analytical eigen-functions corresponding to a steadily propagating crack are used as basis functions for displacements, may move by an arbitrary amount AS in each time-increment At of the numerical time-integration procedure. The moving singular element, within which the crack tip always hgs a fixed location, retains its shape at all times, but the mesh of "regular" (isoparametric) finite elements, surrounding the moving singular element, deforms accordingly. An energy-consistent variational statement was developed in [1,2] as a basis for the above "moving singularity-element" method of fast fracture analysis. It has been demonstrated in [1,2] that the above procedure leads to a direct evaluation of the dynamic K-factor (s), in as much as they are unknown parameters in the assumed basis functions for the singular-element. Solutions to a variety of problems, were obtained by using the above procedure and were discussed in detail in[l,3]. These problems included, among others: (a) constant-velocity, self-similar propagation of a finite central crack in a finite panel [analogous to the well-known Broberg's problem] (b) stress-wave loading of a stationary central crack in a finite panel [analogous to the well-known problems of Baker; Sih et al; and Thau et al.], (c) constantvelocity propagation of a central crack in a panel, wherein, the propagation starts at a finite time after stress-waves from the loaded edge reach the crack, and (d) constant velocity propagation of an edge-crack in a finite panel, whose edges parallel to the crack are subjected to prescribed, time-independent, displacements in a direction normal to the crack-axis analogous to the well-known problems of Nilsson[4]. In Ref.[5], the results of numerical simulation of experimentally measured crack tip vs time history in a rectangular-double-cantilever-beam (RDCB), as reported by Kalthoff et al.[6], were reported. Also (Ref.[5]) the authors' results for the computed dynamic K-factor for RDCB were compared with the experimental (caustics) results of[6], and the independent numerical results of Kobayashi[7], who uses a node-release technique in fast fracture simulation. In this paper, the authors wish first to clarify and comment upon several aspects of the model formulation which appear in their Refs. [1, 2]. Second, the model's accuracy and efficiency are evaluated in terms of less sophisticated models. Finally, the practicality of the special singular element for predicting crack growth for a given crack growth criterion is illustrated. In addressing the second topic, attention is focused on: (a) the effect of using only the stationary-eigen functions (or the well-known Williams' solution) in the moving singularelement for dynamic crack propagation, and (b) the use of isoparametric elements with mid-side nodes shifted so as to yield the appropriate (r -O/2)) singularity[8, 9]. Finally, some recent results are presented which illustrate the facility of the propagation-eigen-function singular element for +Based in part on a presentation made at 2nd International Conference on Numerical Methods in Fracture Mechanics, Swansea, 1980. EFM Vol. 15, No. I-2--M 205 206 T. NISHIOKA et al. predicting crack propagation behavior based on Km vs crack propagation speed as a crack growth criterion. THE VARIATIONAL PRINCIPLE Since the details of the formulation are presented in Refs.[l, 2], only those portions of the formulation necessary to the present discussion will be included here. Further, for simplicity the rather general equations of Refs.[1, 2] will be specialized to the case of Mode I crack growth in bodies subject to traction free crack surfaces, zero body force and with geometry/ applied loading such that the model can make use of symmetry about the crack plane. In Ref. [2], the principal of virtual work proposed as the governing equation for crack propagation during the period [fl, h] in which the crack elongates by amount AE is given as: C O = ] 2 1 2 "2 {(O"/.1" + O'o)~eij + p(Id i + /~il)~//i2} d V Jv 2 fs~2(~ '+ Ti2)6u2dSfa o~jvjlt~lti2dS. (1) The superscripts 1 and 2 refer to quantities at times t~ and t2 respectively. The integrals in order of appearance refer to the volume of the body at time t2, the portion of the surface of the body at time t2 subjected to prescribed tractions and the new crack surface created between times tl and t2. Note that the variational quantities &~ and 8u~ 2 reflect the kinematic constraint at t2 and therefore are arbitrary on AE. Making use of the small strain displacement relation, the symmetry of tr o, and using the divergence theorem, the first term of the volume integral in (1) becomes: fV 2 l 2 -fa (O'2pj2 + O'~jpjl)~ui2 d S f 2 (o-q,s + tr1~,~)au? d V (2) (°'q + ° r q ) 8 ~ i i d V v2 Jv2 2 where ~V2 is the boundary of V2. Noting that (a) AE is part of c~V2 (resulting in the last term of (1) dropping out); (b) that S~ 2 is a part of ~V2; and (c) that 8ui 2 is zero on any portion of ,~V2 that has prescribed displacements, we have: fV 2 "2 1 "'1 2 {tri M -p l l i + trq, i pUi }SUi d V 2 + fs {~2 _ ~ijv 2 + Ti' t~Jv/}~ui 2dS Jar Or2~'2~Ui2 d S = O. (3) Since ~Ui 2 is arbitrary (3) leads to: and or iL j2 = pi~i 2 + fli~i I _ O'#A~ in V2 (4a) o.2p2= ~2+ ~]_o.~,vl on So.2 (4b) 2 2 cr0u i =0 on AE. (4c) It is seen that (4c) is nothing other than the condition that the new crack surface be traction free. Equation (4b) stipulates that traction boundary conditions be satisfied on $~2 and (4a) is a statement of dynamic equilibrium within the body. If one assumes that i = plill then (4a) O ' i j , j reduces to the usual equilibrium equation at t2 (i.e. z _ ..2 2 _ ~rq.s pui ). Then since ~ro. i pzi~ 2 it follows that the state at t3 must also satisfy the usual equilibrium equation and so forth. Therefore, (4a) is equivalent to the standard expression for dynamic equilibrium when the assumption that Moving singularity finite element models 207 1 0-0,~-plii ~ is valid. A similar argument leads to (4b) reducing to the usual condition for satisfaction of traction boundary conditions. In the finite element model formulation, the above assumption is not valid and therefore, the eqns (4) do not reduce to those usually found in finite element model derivations. One reason for 0-~j ~ pti: ~ is that in modeling crack growth, it becomes necessary in the procedure of Refs.[l, 2] to change the mesh configuration at each crack growth time step and to interpolate displacement, velocity and acceleration data at the new node locations. This interpolation will in general result in some disequilibrium in the interpolated solution. If, for example, we assume some disequilibrium at t l such that . . l I _ _ / ~ 2 0 . 2 _ t pui -0 . i i . J [ , then satisfaction of (4a) leads to p ~ 0 . j s . Clearly, the disequilibrium at subsequent steps will be the same form (f) with the sign alternating at each step. Therefore, it would appear that the formulation of Ref.[1, 2] should result in oscillations. Numerical experimentation with the formulation has indeed shown this oscillation to occur. However, when used with the special singular element of Refs.[1,2] the only time it has occurred at discernable levels is in static analyses. In dynamic analyses, it has been found that the inertial forces are generally large enough to make the oscillatory forces negligible. Similar oscillation has been observed when implementing the proposed principle of virtual work with eight-noded isoparametric elements. (Wherein the appropriate crack tip singularity was obtained by shifting midside nodes as suggested by Refs. [8, 9]). It is generally found that the oscillations when using the isoparametric elements are larger than those observed with the special singular element. This is believed to be the result of inherently larger interpolation induced disequilibrium with these less sophisticated elements. A related variational principle for quasi-static crack growth in elastic-plastic bodies was also presented in[l]. To account for effects of history dependent plasticity and finite deformation gradients, an updated Lagrangean rate formulation was used. As a result of the formulation yielding an incremental analysis (as opposed to the computation of total state quantities as in the elasto-dynamic analysis above), the effect of state quantities at tl appearing in the finite element equations for the solution at t2 is fundamentally different. The variational statement, simplified for the case of zero body force, traction free crack surfaces and symmetrical modeling of Mode I crack growth, is given as: [r,itSg 0 Vfs,, ' fazrl, v/ a, dZ =O (5) v, I ' °+S°~g~ld dS+ where ~-lj are the Cauchy stress components of the solution at t, in the current reference configuration (t~, V~, Y~,), ~ii are the incremental displacements in going from the reference state to the final state (t2, V2, E~ +AE) r~ are the applied incremental tractions, g~ =1⁄2(lik.~Uk.j), • L _ I • g,j ~(u;j + tij.~) and S0 are the incremental second Piola-Kirchhoff stress components (or what are also known as Truesdell stress increments). Equation (5), through the use of the divergence theorem, leads to the following Euler-Lagrange equations: l " (rk~Ui.k),~ + S0.j = 0 in V, (6a) (r lk/~i .k+Sis)vi-Ti=O on S~ 1 (6b) . zo)v i =0 on AE. (6c) Equation (6a) is the usual translational equilibrium condition associated with updated-Lagrangean formulations in terms of the second Piola-Kirchhoff stress and does not contain any additional terms such as found in eqn (4a). Equation (6b) is the usual condition for satisfaction of traction boundary conditions and eqn (6c) is the condition that the newly created crack surface be traction free. Further study of the derivation of eqns (5) and (6) shows that there is nothing inherent in the formulation to account for (or correct) the error from interpolation of quantities from the finite element mesh at t~ with crack length E~ to the finite element mesh at t~ with crack length YI+AE. However, the accumulated error at increment P(~p) can be measured by: ~p = ro ~gij d V T/eSa; dS. p 208 T. NISHIOKA et aL which is a check on the equilibrium at increment p. Since this rate formulation does not result in terms which lead to oscillation of the solution it seems that formulation of the linear-elastic dynamic problem in terms of an incremental model similar to that for the elastic-plastic large deformation problem would eliminate the trouble with oscillations while at the same time retain the crack growth modeling features. Results of analyses which are presented later in this paper are largely based on the formulation as originally proposed in Refs.[1,2] since it appears that no significant change would result from a reformulation. The one exception to this are the results obtained through the use of isoparametric elements with midside nodes shifted so as give a singularity at the crack tip. In those analyses, the usual statement of virtual work is applied: : { O ' i j ~ i i " + flU i t~ll i } d V Ti2~ui 2 d S . 2 cr 2 Traction free crack surfaces are approximated by letting nodal forces on the crack surfaces be zero. SINGULAR ELEMENT FOR DYNAMIC CRACK PROPAGATION A singular crack tip element was also developed in[l, 2] and used in conjunction with the formulation (1) for the analysis of Mode I dynamic crack propagation in linear elastic two dimensional bodies. This singular element uses an arbitrary number of the displacement eigen-functions which come from the solution of a crack in an infinite body. For dynamic crack propagation, these eigen-functions are taken as those for the corresponding steady-state dynamically propagating crack in an infinite body. Equations (7) give the form for the assumed displacement, velocity and acceleration within the singular element: u'(~:, x2, t) = U(~:, x2, v)13(t) (7a) ti s = Ufl v(U), J] (Tb) iP : U/~ 2v(.U),efl + v2(U),~d] (7c) where U is the matrix of eigen-functions (plus appropriate rigid body modes), fl is the vector of undetermined coefficients, v is the crack propagation speed and (sclX2) are the coordinates relative to the moving crack tip (~ = x~ vt) . It should be noted that (7b) and (7c) are obtained from (7a) through differentiation with respect to time with the assumption that v is not a function of time. The following comments should remove any incorrect notions that this in any way limits the use of the element to constant speed crack propagation. First, it has been shown [12] that the near-tip fields are the same for steady-state and transient crack propagation. Therefore, provided v at each time step reflects the current speed, there is no question that the eigen-functions for the element are correct and that the coefficient of the singular eigenfunction (/31) is indeed the Mode I stress intensity factor. A second consideration is that the associated stress eigen-functions do not satisfy the stress equilibrium-equation for non-steadystate (as viewed by an observer moving with the crack tip). ['821li _ (~2U i 2 (92//i ~ = T)o (8) but instead, satisfy the corresponding steady equation: 2 02ui ~o,J = pv ' ~ . (9) While it would be preferred that eqn (8) be satisfied exactly, the displacement finite element method does not require this. The stress equilibrium of (8) is satisfied in the usual approximate sense associated with the finite element method. One common difficulty which arises when using more than one element type in a model is Moving singularity finite element models 209 the lack of compatibility at the interface of the dissimilar elements. In Refs.[1,2] a method is explained for maintaining compatibility at the boundaries of the singular element which are shared by eight-noded isoparametric elements. This method involves selecting/I such that the following three error functionals are minimized: I~=fp (u~-uR)2dp; I2=fp 0PfR)2dp ; I3=f (ii '-iiR)2dp (10)